from scipy.special import gdtrix, ndtri
# import matplotlib.pylab as plt
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from sklearn import metrics
import time
from sqlalchemy import create_engine

# 解决matplotlib显示中文乱码问题
plt.rcParams['font.family'] = ['sans-serif']
plt.rcParams['font.sans-serif'] = ['FangSong_GB2312']  # windows


# plt.rcParams['font.sans-serif'] = ['Arial Unicode MS']  # Mac


def xtrans(plist):
    """
        x坐标轴转换成海森机率格纸格式
        :param plist: plist传入的是一个p值列表数据，plist必须是概率值，一定是大于1的数字，如6%，取6，返回的xplist也是一个列表数据
        :return xplist: 符合海森机率格式的x坐标轴
        """
    xzero = ndtri(0.0001)
    xplist = []
    for i in range(len(plist)):
        xplisti = ndtri(plist[i] / 100)
        xplisti -= xzero
        xplist.append(xplisti)
    return xplist


def ljxs(plist, Cs):
    """
        求离均系数fp
        :param plist: plist传入的是一个p值列表数据，plist必须是概率值，一定是大于1的数字，如6%，取6，返回的xplist也是一个列表数据
        :param Cs:
        :return:
        """
    aa = 4 / Cs ** 2  # b值，反映的是形状（shape），在水文计算中为α（alpha），不要与scipy中的官方文档中的a值混淆
    tp = []  # 为求离均系数fp，tp是过渡
    fp = []
    for i in range(len(plist)):  # 离均系数fp
        tpi = round(gdtrix(1, aa, 1 - plist[i] / 100), 3)  # 官方文档中的a值为scale，本次采用了标准伽马分布，所以取a=1
        fpi = round(Cs * tpi / 2 - 2 / Cs, 3)
        tp.append(tpi)
        fp.append(fpi)
    return fp  # 返回的fp也是一个列表数据


# 坐标轴的绘制

def draw(file, imagename="", N=2, rain_dim=""):
    # file = r'/home/matrix/workspace/StudyNotes/暴雨洪涝研究/气候中心-P3/inputFile/福州.txt'
    global Q, df, sum_nan
    try:
        if file[-3:] == 'txt':
            df = pd.read_table(file)
        elif file[-3:] == 'xls' or file[-3:] == 'lsx':
            df = pd.read_excel(file)
        df.columns = ['col1', 'col2']
        df = df.sort_values(by='col2', ascending=False, axis=0)
        sum_nan = df.isnull().sum()  # 统计缺失值总数
        mean_pre = df.iloc[:, 1].mean()
        df.fillna(value=mean_pre, inplace=True)
        Q = df.iloc[:, 1].astype("int")
    except Exception as e:
        print(e)

    pstandard = [0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1.0, 2.0, 3.0, 5.0, 10.0, 20.0, 25.0,
                 30.0, 40.0, 50.0, 60.0, 70.0, 75.0, 80.0, 90.0, 95.0, 97.0, 99.0, 99.5, 99.9]
    # pstandard = np.linspace(0.01, 100, len(Q)).round(2)
    x_axis = [str(i) + '%' for i in pstandard]
    x = xtrans(pstandard)
    y_axis = np.linspace(0, round(max(Q) * 1.2, 2), 10).round(1)
    # print(x_axis, y_axis, sep='\n')
    # print(x)
    # 建立海森图纸张，需要知道Q的最大值来确定横坐标的上限
    fig = plt.figure(figsize=(12, 6))
    ax1 = fig.add_subplot(111)
    for i in range(len(x)):
        plt.vlines(x[i], 0, round(max(Q) * 1.2, 2), 'g', '--')
    for j in range(len(y_axis)):
        plt.hlines(y_axis[j], 0, max(y_axis), 'g', '--')
    plt.xticks(x, x_axis, color='blue', rotation=90)
    plt.yticks(y_axis, y_axis, color='blue', rotation=0)
    plt.ylim(0, round(max(Q) * 1.2, 2))
    plt.xlim(0, round(max(x) + 0.1, 2))
    plt.tick_params(labelsize=10)
    ax1.set_title("{}".format(imagename))
    ax1.set_xlabel("概率设定值P%")
    ax1.set_ylabel("降水量{}".format(rain_dim))

    # 绘制原始数据的散点图
    total = len(Q)  # 样本总数
    Ex = np.mean(Q)  # 均值
    qcv = np.std(Q) * np.sqrt(len(Q) / (len(Q) - 1)) / Ex  # 无偏估计的Cv值,标准差除均值
    print(np.std(Q))
    # Q.sort(reverse=True)
    # print("降水量Q（0.1mm）:\n {}".format(Q))
    # Q的经验频率 p=m/(n+1)，这里乘以100的原因是pstandard取值为（0.01~99.9）
    pq = pd.Series([(i + 1) * 100 / (len(Q) + 1) for i in range(len(Q))])
    # print("Q的经验频率 p=m/(n+1):\n{}".format(pq.astype("int")))
    Pq = 100 / pq  # 重现期 P = 1/p，为了画图的方便，之前pstandard乘以100，故此处用100/pq
    pqx = xtrans(pq)
    plt.scatter(pqx, Q)
    # 绘制预测的概率曲线，需要知道Cs,Cv,Ex（均值），一般通过倍比n=Cs/Cv来计算Cs=Cv*n
    n = N  # 是自定义值
    Cs = n * qcv
    fp_real = ljxs(np.linspace(0.01, 99.9, len(Q)).round(2), Cs)
    fp = ljxs(pstandard, Cs)  # 离均系数（画图用）
    qpredict = [Ex * (qcv * m + 1) for m in fp]  # 设计值xp
    qpredict_real = [Ex * (qcv * m + 1) for m in fp_real]  # len(Q)设计值xp
    RMSE = np.sqrt(metrics.mean_squared_error(pqx, qpredict_real))  # 绝对均方误差
    # print("fp:\n {}\n qpredict=:\n{}\n Cs = {} \n qcv = :\n{}\n pqx=:\n{}".format(fp, qpredict, Cs, qcv, pqx))
    plt.plot(x, qpredict, 'r-', lw=3)
    plt.text(6.5, max(Q) * 1.1, r'mean_value={:.2f},Cs={:.2f}, Cv={:.2f}'.format(Ex, Cs, qcv), size=15,
             family="Monospace", color="black", weight="normal", verticalalignment='top',
             horizontalalignment='right', bbox=dict(facecolor="w", alpha=1))
    # plt.show()
    tempimage = "temp.png"
    plt.savefig(tempimage)
    # return saveimage
    return total, sum_nan, Ex, Cs, qcv, fp_real, qpredict_real, Pq, RMSE, tempimage
